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Abstract 

We present in this paper detailed numerical Vlasov simulations of the Hamil- 
tonian Mean-Field model. This model is used as a representative of the class 
of systems under long-range interactions. We check existing results on the sta- 
bility of the homogeneous situation and analyze numerical properties of the 
semi-Lagrangian time-split algorithm for solving the Vlasov equation. We also 
detail limitations due to finite resolution of the method. 



1. Introduction 

Numerical resolution of the Vlasov equation plays a crucial role in the field 
of plasma physics, as a theoretical investigation tool or to study realistic plasma 
devices. Development of Vlasov codes is still an active research topic [Ij, pre- 
senting technical aspects (parallelization, memory usage limitations, performing 
full 6D simulations 3|) as well as theoretical ones (numerical dissipation, "non- 
Vlasov" effects (J, Hi). Recently, the Vlasov framework has been used to study 
toy models of particles under long-range interactions 0, S S Hi ■ These models 
comport an all to all coupling and model more complex physical situations such 
as Coulomb or gravitational interactions. A peculiarly slow relaxation to ther- 
modynamical equilibrium has been discovered in these systems, and understood 
with the help of the Vlasov equation. 

A model displaying many phenomena related to long-range interactions is 
the Hamiltonian Mean-Field model (HMF) ^ . It displays an equilibrium phase- 
transition from an homogeneous to a non-homogeneous state, characterized by a 
mean-field order parameter. Depending on the initial condition, the HMF expe- 
riences an out-of-equilibrium phase transition. A theory developed by Lynden- 
Bell for galactic collisionless dynamics predicts the coarse-grained asymptotic 
evolution of the particles' distribution. It has been succesfuUy applied to the 
HMF, predicting the out-of-equilibrium phase diagram [1, 

Only one study related to the HMF makes use of a Vlasov simulation 7] . A 
good knowledge of the numerical procedure, applied to models with long-range 
interactions, is however necessary if one aims at widening the toolset used to 
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study their dynamical properties, especially the concept of the so-called Quasi- 
Stationary States that depend on the number of particles considered. The other 
way round, the use of simple models in Vlasov simulation may offer a better 
analysis of the tool itself, serving the needs of plasma physics research. 

We present in this article detailed numerical simulations of the Vlasov equa- 
tion associated to the HMF model. We check existing results obtained via 
N-body simulations and analyze the numerical aspects of the simulations. 

Section O presents the HMF, section [3] presents the algorithm used. We 
present the results of the computations in section [J] and conclude in section [SJ 



2. Modeled system 

The Hamiltonian Mean-Field (HMF), introduced in Ref. [§], consists of N 
rotators interacting with a cosine attractive potential. It provides a simplified 
model of self-gravitating sheets, its potential corresponding to the first Fourier 
mode of the latter model. The Hamiltonian is : 

N 2 1 ^ 

^ = Ey + ^ E (i-cos(0.-0,)), (1) 



= 1 ij = l 



we define the vector : 

M : 



M = |M| is called the magnetization and quantifies the degree of inhomogeneity 
of the particles; it also gives a simplified expression for the interaction part of 
the Hamiltonian, V = |(1 — M^). We shall make use of M to follow the time 
evolution of the system. 

In the mean-field limit N ^ oo, the HMF is described by a Vlasov equation : 

5/ df dv[f] df _ 

dt ^de de dp 

V[fm = l-AUf]cos0-My[f]sine 

M,[f] - J dOdpf cose, 

M,[f] = [dOdpf sine (3) 



where / = f{e,p;t) is the one-particle distribution function, V is the potential, 
depending self-consistently of /. The energy is conserved, and takes the form 

Uit)[f] ^ Jdedp fie,p;t) (^^ + i (1 - cose- My[f] sine)^ . (4) 

In the following, we make use of homogeneous initial conditions, f{e,p) = 
f{p) and study their consequent dynamics. An homogeneous state is by defini- 
tion stationary because and My (hence the force-field) are zero, but may be 
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stable or unstable. A stability criterion has been devised by Yamaguchi et al [5| 
(see also Chavanis jH). Waterbag and Gaussian distributions display an energy 
treshold below which an unstability occurs, giving rise to a magnetized state; 
N-body simulations were performed in Q to check these predictions. The mag- 
netized state displays oscillations of M around a value that is found to be below 
thermodynamical equilibrium. Finite N effects eventually allow the system to 
reach equilibrium. 

Lynden-Bell has devised a statistical theory to predict the outcome of the 
coUisionless galactic dynamics. It has been successfully applied to the HMF, 
in the case of waterbag initial conditions and compared to Vlasov simu- 
lations 0. It predicts a / that is different from equilibrium; this solution is 
interpreted as an intermediary stage of the N-body dynamics, hence the term 
Quasi-Stationary State, whose lifetime diverges with N . 



3. Algorithm 

We use the algorithm described in Ref. reviewed with cubic spline 

interpolation [11]. It is based on the characteristics method for solving a partial 
differential equation : the value of f{0,p) at time t + At is taken from / at time 
t, at the foot of the characteristic curve: 

r+\e,p) = ri0 ~ At{p+^F*{0)At),p~ AtF*i0)) (5) 

where 6 = 6 — pAt/2 and F* is the force field computed at half a time-step. We 
use the notation f^{0,p) = f{9,p;s At). 

A practical implementation to compute Eq. ([5]) on a numerical mesh is the 
following : 

1. Advection in the f*{e,,Pm) = ^"{0^ - Pm.At/2,pm) 
0-direction, 1/2 time-step 

2. Computation of the force field 
for /* 

3. Advection in the f**{0^,Pm) = f*{0,,p„, ~ F*{9,)At) 
p-direction, 1 time-step 

4. Advection in the f'^H^i^Pm) = - PmAt/2,p„) 
^-direction, 1/2 time-step 

Each of these steps is performed for the whole grid before going to the next 
one. A sketch of these steps is given in Fig. [TJ This method, involving a fixed 
mesh and trajectories along the characteristics backwards in time is called semi- 
Lagrangian. It is used in fluid dynamics, for weather forecasts for instance (see 
Ref. [i3). 

The time-split approach is valid at the second order for the Vlasov equation 
[loj . An interpolation scheme must be provided in order to compute .P(^, p) 
outside of the storage mesh; we will use cubic spline interpolation [l3|. The 
interpolation needs only to be 1-dimensional, whether in the 0-space or in the 
p-space, simplifying the numerical procedure. Cubic splines are commonly used 
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in plasma simulations, but other advection algorithms exists that can improve 
conservation or monotonicity properties, depending on the situation, see Refs. 



a) b) c) 



Figure 1: The semi-Lagrangian algorithm: a) illustrates the computation of the value of / at 
the grid points indicated by the •'s from the feet of the characteristic curves (the start of the 
arrow), b) and c) represent steps 1 and 3 of the time-split algorithm. 



Altough the Vlasov equation conserves theoretically every invariant of the 
form Cslf] = J dO dp s{f{d,p)) (the so-called Casimirs), the numerical scheme 
shows deviations. We define the Li(t) norms : 



L,(t)= / d0dpifi9,p)y 



(6) 



Li is the mass of the particles (it is normalized to 1); we will use L2 to charac- 
terize numerical dissipation. 

The computations were performed in Fortran 90, and the data read with the 
help of numpy, hSpy and matplotlib 1^, 17, T^ . 



4. Results 

We discuss in this section the computations related to the properties of the 
simulation for an initial gaussian homogeneous distribution, then we check the 
stability criterion of Ref. We also compare the runs in the spirit of Refs. 

4-1. Properties oj the numerical solution 

We run simulations starting with a gaussian profile, slighty perturbed : 

f{e,p)^\[^e-^P'/^\ + es\ne) (7) 

at an energy oiU = 0.51, e — 10^^. (3 is computed from U , for an homogeneous 
state : /3 = | u-1/2 • At J7 = 0.51, the gaussian profile is unstable; M grows, 
then oscillates around a mean value to which it eventually relaxes. 
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The following parameters are used : Ai = 0.1, size of the box in the p- 
direction [—4.5 : 4.5]. The resolution is varied from Ng = Np — 64 (G64) to 
Ne = Np = 512 (G512); we call these runs G64, G128, G256 and G512. 

We detail in Table [T] the conservation properties of the algorithm. The 
conservation of U could be improved by decreasing At, implying more compu- 
tational time. This is the sole constraint on At as the semi-Lagrangian method 
has no Courant condition. 



run 


Li(0) 


max ^(.Q^ 


G64 


1.6 10-4 


3.2 10~4 


G128 


2.1 10-5 


2.0 10-4 


G256 


7.2 10-6 


2.0 10-4 


G512 


2.8 10-6 


2.0 10-4 



Table 1: Conservation properties at different resolutions. 
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Figure 2: Evolution of M for the gaussian ini- Figure 3: Evolution of L2, same runs as Fig. [2] 
tial condition, Eq. lO (f/ = 0.51,e = 10"*). 
A magnification of the same curve is presented 
in the inset. 

M displays a similar behaviour for all values of the resolution until t sa 60, 
except for the lowest resolution run G64; M has been checked to correspond to 
N-body dynamics with N = 10^ for short times (data not shown). M grows 
up to a saturation value, identical for all runs, then oscillations take place and 
eventually damp to a stable value. 

The damping of the oscillations is strongest in G64, they disappear at i « 70. 
This damping is present in the other runs as well, a phenomenon that diminishes 
when the resolution is increased. This fact is to be put in correlation with the 
evolution of L2 (Fig. [3]) : The decrease of L2 means that the small scales are 
not well described anymore by the numerical procedure. This decrease is the 
strongest in G64 where L2 doesn't reach a plateau and diminishes continuously. 
The shape of /, a vortex rotating in phase space, is giving the oscillations of M. 
When this structure is smoothed out by the lost of the small scales, see Fig. El 
M tends to a constant value. Run G512 displays a filamentary structure and 
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oscillations in M up to the end of the simulation. 

Asymptotically, all runs tend to a similar value of M, which we understand 
with the standard deviation of M, ctm, and the mean value across simulations, 
M. <tm/M = 1.1%, or 0.7% if we leave G64 aside. 

4-.2. Stability of the initial condition 

The stability of homogeneous initial condition has been studied in . We 
perform a check for the waterbag and gaussian initial profiles. The waterbag is 
the following profile, for a given U: 

= { "^^^ (8) 

Starting from both side of the critical energy {U* — 7/12 for the waterbag 
and U* = 3/4 for the gaussian), we observe in Figs. 0] and [5] the change in 
stability. For values of U above U* , M stays at the same order of magnitude or 
decreases. 




Figure 4: M (t) for initial gaussian profiles. M 
grows if U is below U* = 3/4. 



Figure 5: Af(t) for initial waterbag profiles. M 
grows if U is below U* = 7/12 ?5! 0.58. The 
slope before saturation of the three first runs is 
used to compute the exponential growth factor 



In the case of the waterbag, we compared the exponential growth rate de- 
fined in Ref. [3] with the one measured from the runs of Fig. [5] and found a 
very good agreement. The exponential growth rate depends on the energy U : 
A = y/ 6 (U* — U) . The theoretical and numerically computed values are given 
in couples (A^jj,Anum) for U = 0.52,0.55 and 0.57 respectively : (0.61,0.61), 
(0.43,0.43) and (0.32,0.31). 

4-3. Effect of the resolution 

The results of section 14.11 considered different resolutions of the same prob- 
lem; f{0,p) was inspected visually, but most of the analysis focused on macro- 
scopic quantities : Li, L2, U and M. Galeotti, Califano and Mangeney [1,0] 
compared cuts in /, i.e. f{9,p*) where p* is fixed or f{9*,p) where 9* is fixed, 
after smoothing of all runs to the lowest resolution. This method explores with 
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more precision the structure of /. For runs G64 to G512, we display f{9, 0) in 
Figs. [6]and[7]at two times : t — 6A when runs G128 to G512 still have a similar 
M, and at t = 400 in order to discuss the asymptotic evolution. 




-3 -2-10 1 2 3 -3 -2 -10 12 3 
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Figure 6: Cut in the function / at time t = 64. Figure 7: Cut in the function / at time t = 
Same runs as Fig. [2]The locations of the peaks 400. Same runs as Fig. [2] The location of 
is similar in all runs. the maximum is the same in this asymptotic 

state. 

Ki t — 64, altough L2 values are different, which may imply a different 
evolution of /, we find peaks in f{0,0) at the same locations. This indicates 
that / behaves similarly. Indeed, the energy spectrum of the HMF has only 
one component and the force on the particles (or phase space elements) only 
depends on M. Thus, as long as M is identical for the runs, so is the force term. 

At t = 400, G64 is completely smooth, and the other runs still show a small 
filamentary structure. The location of the maximum is the same, a fact that 
contrasts with the results of Ref. for a Vlasov plasma, and that indicates a 
similar asymptotic state. 

We gain by these observations, for this model, the fact that one observable, 
M, not only describes the macroscopic evolution, but by comparison between 
runs also monitors the structure of /. Development of Vlasov codes may benefit 
from this simplicity while comparing different algorithm, a process often used 
in numerical analysis (for Vlasov plasmas, see for instance Ref. [3]). 

5. Conclusions 

We presented in this paper a detailed study of Vlasov numerical simulations 
for the Hamiltonian Mean-Field model (HMF), that have been compared to 
previous work on the HMF and to Vlasov results in plasma physics. We gave 
evidence that the mean-field order parameter M of the HMF is sufficient to 
compare different runs in terms of phase space evolution as well as dissipative 
(i.e. damping) properties. 

Illustration was given of the "non- Vlasov" limit of the numerical simulation 
which we need to keep in mind to analyze simulations results carefully, as stated 
recently by other authors d, 0, 0] ■ 

Finally, we hope that cross research between long-range systems obeying a 
Vlasov equation and Vlasov plasmas may benefit to both fields by providing the 
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Figure 8; Phase space contour lines, a) and b) are from run G64 at times 32 and 64. c) and d) 
are from run G512 at the same times. While a) still presents the correct shape as compared 
to c), b) has lost most of the spiral shape present in d). 

first with a new tool they can apprehend and the latter with new insight on the 
numerical method. 
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